Experimentally switching between liquid and solid phases?

Looking at Grochola's free energy path, I thought it was a nice computational technique.

The March 2004 issue of Physics Today has an article on controlling atoms and ions in optical traps (it's actually about using these for quantum information processing). On page 41, it describes putting the atoms in a optical lattice. The depth of the wells and the hopping matrix elements can be controlled by varying the intensity of the lasers creating the lattice.

This seems like one could experimentally switch between solid and liquid states without going through a phase transition.

There are some differences - I don't know if they are important or not. First, the experimental setup is a quantum system, Grochola's paper describes a classical system. Second, the confining potential is different. In the experimental setup, atoms are confined to any lattice site, but are not tethered to a particular lattice site.

The Physics Today article mentions a Mott insulator transition when the confining potential becomes strong enough. This transition seems due to the superfluid (quantum) nature of the liquid state. If the experiment were performed on a non-superfluid system, this transition would not happen ??

Update (3/16/2004): I forgot the final stage in the path - turning off the confining potential to get a stable solid state. I'm guessing the lattice spacing of the optical lattice is at least the laser wavelength or greater, so it's much larger than the lattice spacing of a normal solid. But some atoms can have their scattering length tuned - could it be made large enough to be comparable to the optical lattice spacing? In that case the atoms would behave as hard spheres (most importantly, there would likely not be enough attraction to keep the solid together), so an external potential would be required. Most of these setups seem to generate a harmonic confining potential (a box-like potential seems much harder to generate??), so the question is - what is the behavior of hard spheres in a harmonic potential?

For this system, there are two parameters - the number of particles and the strength of the external potential. I suspect there would be three regimes - one all liquid, one with a solid center surrounded by a liquid exterior, and one would be all solid. To take the thermodynamic (large N) limit, it seems that one would need to increase the particle number and decrease the external potential strength simultaneously to keep the central density or pressure constant (otherwise the pressure/density at the center would increase with added particles). Or maybe such a limit doesn't exist?

The Way of Antisymmetry

The traditional QMC wavefunction is the product of a determinant of single particle orbitals and a Jastrow factor (which contains the electron-electron correlation). The determinant ensures the required spatial antisymmetry of the wavefunction, and is computable in O(N³) time (where N is the number of electrons). But it is not the most general form possible.

The most general form is obtained by explicitly antisymmetrizing: summing the values of a function with all possible permutations of the electron coordinates (and multiplying by the sign appropriate for the permutation). For example, with 2 particles and a function f(r₁,r₂), the antisymmetrized form is f(r₁,r₂) - f(r₂,r₁).

Unfortunately, the number of terms is N!. (The determinant form is obtained by applying this procedure when the function has the special form of the product of functions that each depend only on a single electron coordinate - these functions are the single particle orbitals. In the two particle example, f = f₁(r₁)f₂(r₂). The antisymmetrized form is f₁(r₁)f₂(r₂) - f₁(r₂)f₂(r₁), which is the determinant of a matrix.)

Dario Bressanini, Massimo Mella, Gabriele Morosi, and Luca Bertini have several papers using this form of the wavefunction for small molecules, and they get very good results. The papers can be found on Massimo Mella's publication list (edit 12/3/2005 - removed nonworking link). They are (this may not be a complete list):

• Many-electron correlated exponential wavefunctions. A Quantum Monte Carlo application to H2 and He2+
  
  Chemical Physics Letters 240, 566 (1995)
  
  Paper available at this Citeseer page
  
  
• Nonadiabatic wavefunctions as linear expansions of correlated exponentials. A Quantum Monte Carlo application to H2+ and Ps2
  
  Chem. Phys. Lett. 272, 370 (1997)
  
  Paper available at this Citeseer page
  
  
• Linear expansions of correlated functions: a variational Monte Carlo case study
  
  Int. J. Quantum. Chem. 74, 23 (1999)
  
  Paper available at this Citeseer page
  
  
• Explicitly correlated trial wave functions in quantum Monte Carlo calculations of excited state of Be and Be-
  
  J. Phys. B, (2001)
  

Not all N! terms will give a significant contribution to the result. Is it possible to use techniques similar to O(N) approaches to reduce the number of terms? Is it possible to sample terms from the sum? These might have trouble near the nodes, since the wavefunction vanishes due to the cancelation of all the terms.

DMC of porphyrin

An all electron calculation of a large molecule (large by QMC standards, anyway):
Quantum Monte Carlo for electronic excitations of free-base porphyrin by Alan Aspuru-Guzik, Ouafae El Akramine, Jeffrey Grossman, and William Lester. The porphyrin they used has the formula C₂₀ N₄ H₁₄, so 168 electrons.

Solid vs. Liquid, The Free Energy Battle

Determining the free energy difference between the solid and liquid states is essential for accurately mapping the phase diagram of a substance. Typically, one determines the free energies of the solid and liquid states through separate thermodynamic integrations, where each integration has one end point connected to a state with an analytically known free energy. In the solid state, the Einstein crystal (springs attaching particles to lattice sites) is often used as a starting point. In the liquid state, the low density gas can be used as an end point. Other more complicated scalings of the potential are also possible for the integration path.

A direct integration path between the solid and liquid states would be useful.
Gregory Grochola has published one such path in "Constrained fluid lambda-integration: Constructing a reversible thermodynamic path between the solid and liquid state"(JCP 120, 2122).

The basic path is to start with the liquid state and turn off the intermolecular potential. Then turn on a potential constraining particles to lattice sites (He used a Gaussian potential). The last step is to turn off the lattice potential while turning the intermolecular potential back on.

He mentions at the end that this method is not very computationally efficient. I would expect using Bennett's method for free energy differences would cut down on the number of intermediate states needed.

Questions

1. Compare and contrast this path with the phase switching Monte Carlo technique (there's a reference in the article)
   
2. Thermodynamic integration requires there be no phase transitions on the path. Is it immediately obvious that this path has no transitions? Under what conditions, if any, would a phase transition occur on the path?
   

Hard Spheres on the Move

Cluster algorithms increase efficiency for systems on a lattice by flipping large numbers of connected spins at once. In "Cluster Monte Carlo algorithms", Werner Krauth summarizes work on how such algorithms can be used in a continuum hard sphere simulation.

The basic idea is to rotate particles by &pi around a randomly chosen axis. (The key observation is that this transformation is its own inverse). A practical algorithm (the pocket algorithm) chooses a collection of particles that can be rotated without any hard spheres overlapping. It starts with a single particle and rotates it. If it overlaps with any other particles, they are added to the collection. Then the new particles are tested in the same way, and the collection is built up of all the particles that need to be rotated to prevent any overlaps.

At high densities, this algorithm will fail because the pocket cluster will contain all (or most) of the particles in the system, and simply rotate the entire configuration. Despite this limitation, there are systems where these moves are useful (see the paper for some examples).

On another subject, suppose you wanted to keep the trial box size fixed, and see how the dynamics varied with density (which would be useful if one were studying the actual MC dynamics) At higher densities, the efficiency would get quite bad. Hirosi Watanabe, et al tackle this problem in A rejection-free Monte Carlo method for the hard-disk system. They use geometry rather than random sampling to determine the allowable move area (and acceptance probability). To actually choose a point, they use a rejection method (it's not clear why this method is called "rejection-free", then. Maybe because the rejections are not essential for maintaining detailed balance, and rejected moves don't need to be included in the average.). At the step of choosing a new point, the area they choose from is greatly reduced from the original trial area.

Reducing the Variance

Assaraf and Caffarel have published a couple of papers (JCP 119, 10536 and JCP 113, 4028) about reducing the variance for general estimators in QMC. The VMC energy has the well known zero variance property - the variance goes to zero as the trial wavefunction approaches an eigenstate. Through a "renormalization" , the authors show how to extend this property to quantities other than the energy.

The basic approach seems to be to add a term (with zero average) to the bare estimator. Then that added term can be optimized to reduce the variance. Their derivation starts with the Hellmann-Feynman theorem, but this seems confusing to me.

Why not start with the new estimator, O tilde = O + B. ("O tilde" and "O" are notations from the paper. "B" is my own notation for representing all the added terms) Then we demand < B > = 0, and use a Lagrange multiplier to enforce it while optimizing the variance (<(O+B)^2> - < O+B>^2)?

Their paper gives a particular form for B, that ensures the < B>=0 constraint is met (or at least that < B > vanishes as the trial wavefunction approaches the ground state).

Can this approach be used to reduce the variance in classical Monte Carlo as well? In principle, B could optimized in that case as well - the hard part is finding a special form as in Quantum Monte Carlo.

Added thoughts (12/14/2003): The primary application of this method is to forces, which are plagued by an infinite variance in QMC. This looks a like good way to improve computations of the force.

The paper looks at average values and variances as functions (or at least the polynomial order) of delta psi (the diference between the exact ground state and the trial wavefunction). For general (non-energy) estimators, the bias (or error) in the average value is linear in delta psi, and the variance is constant. Using their procedure, the error is quadratic and the variance is linear. It appears the added terms (B) cancel out the leading error terms in the estimator. The question is then, why stop at the leading order? Can a similar process by applied to cancel out the quadratic term in the average ??

Free Energy Chronicles

Free energy calculations are annoyingly difficult. Computing averages with Metropolis is relatively easy,
but the free energy is hard to get. The usual method is thermodynamic integration over a path from a state of known free energy to the desired state. (Maybe it's the ensemble averages that are anomalously easy, making other quantities seem difficult by comparison) In addition, some other aspects about free energy calculations are bothering me as well.

The typical example for explaining Monte Carlo integration is to compute the area of an irregularly shaped field or pond by throwing rocks at it. But then computing the area (volume, free energy, whatever) is actually quite difficult in practice. Of course in this example, one has to compare the rocks falling in desired area to those falling in an enclosing area that can be computed analytically, so it's still computing the difference between two states.

Another issue that bothers my intuition is that if I were to wander around in a field I would have some idea how big the field is based on my wanderings. In other words, why can't we get the area from properties of the random walk. This could simply be pushing the analogy too far. A person wandering is not undergoing a truly random walk. They would likely walk in single direction until hitting an edge and then change direction. So imagine riding an out-of-control Segway around a field - would you still get a sense of the size of the field.

Additionally, a person is able to see the edges. So imagine being in a corn field (I am in central Illinois, after all - others can imagine a forest, with very soft trees) on an out-of-control Segway. Would it still be possible to get an idea of how big the field is based on the path taken and how often you bump into the edge?

So I'd like to look into some of the properties of random walks - how much information can be extracted. Are there theoretical limits? Can we change the random walk or come up with different estimators that would let us compute the area?

In order to investigate concrete but simple systems, I'm thinking of a series of systems to make the problem increasingly realistic for physical systems. The simplest and easiest to think about the area of an irregular region. More realistic, the area of the region is much less that the area of the enclosing square, and the region is long and snaky. Then go to higher dimensions, and eventually to a system of hard spheres. To be more generally applicable, the next step would be a general potential (think of a hilly field and we want to know the volume of soil is in the field)